##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
suppressWarnings(library(edgeR, quietly = T))

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/RSEM/CombineCounts.FilterLowExpression-MergeMutiSample.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_tpm <- data.frame(fread(rsem_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

##########################################################################################

class_type <- c( "Normal" , "IM" , "IGC" , "DGC")

##########################################################################################

dat_tpm_all <- dat_tpm[,-2]

##########################################################################################
## 构建分组信息
condition <- factor(sapply(strsplit(colnames(dat_tpm_all)[-1] , "_") , "[" , 2))    
#individual <- factor(sapply(strsplit(colnames(dat_tpm_all)[-1] , "_") , "[" , 1))
#coldata <- data.frame(row.names = colnames(dat_tpm_all)[-1] , condition , individual)

## counts变为整数
readCount <- round(dat_tpm_all[,-1])
rownames(readCount) <- dat_tpm_all$gene_id

##########################################################################################
## https://rpubs.com/LiYumei/806213
conditions <- condition

## Count matrix preprocessing using edgeR package
y <- DGEList(counts = readCount , group = conditions)

## 去除TPM低表达的基因
## Remove rows consistently have zero or very low counts
#keep <- filterByExpr(y)
#y <- y[keep,keep.lib.sizes=FALSE]

## Perform TMM normalization and transfer to CPM (Counts Per Million)
## 统一标化
y <- calcNormFactors(y,method="TMM")
count_norm <- cpm(y)
count_norm <- as.data.frame(count_norm)

##########################################################################################
## 标准化的矩阵,用于后续样本间的表达比较
## 后续除软件要求计算评分使用TPM，否则均使用该矩阵
normalized_counts <- count_norm
normalized_counts <- data.frame(normalized_counts)
normalized_counts$gene_id <- rownames(normalized_counts)
image_name <- paste0( out_path , "/CombineCounts.FilterLowExpression-MergeMutiSample.TMM.tsv" )
write.table( normalized_counts , image_name , row.names = F , sep = "\t" , quote = F )


##########################################################################################
## 用于加权神经共表达网络
## https://www.biostars.org/p/280650/
dat_forwgcna <- log2(count_norm+1)
dat_forwgcna_out <- data.frame(dat_forwgcna)
dat_forwgcna_out$gene_id <- rownames(dat_forwgcna_out)

image_name <- paste0( out_path , "/CombineCounts.FilterLowExpression-MergeMutiSample.TMM_log2.tsv" )
write.table( dat_forwgcna_out , image_name , row.names = F , sep = "\t" , quote = F )


##########################################################################################
## 多组，两两做差异表达
result_diff <- c()

for( k in 1:(length(class_type)-1) ){
    class1 <- class_type[k]
    print(class1)

    for( j in (k+1):length(class_type) ){
        class2 <- class_type[j]
        print(class2)

        ## Run the Wilcoxon rank-sum test for each gene
        use_sample <- grep( paste0(class1 , "|" , class2) , colnames(count_norm) , value = T)
        count_norm_use <- count_norm[,use_sample]
        conditions_use <- grep( paste0(class1 , "|" , class2) , conditions , value = T)

        pvalues <- sapply(1:nrow(count_norm_use),function(i){
            data<-cbind.data.frame(gene=as.numeric(t(count_norm_use[i,])),conditions_use)
            p <- wilcox.test(gene~conditions_use, data)$p.value
            return(p)
        })
        fdr=p.adjust(pvalues,method = "fdr")

        ## Calculate the fold-change for each gene
        conditionsLevel <- c(class1 , class2)
        dataCon1 <- count_norm_use[,c(which(conditions_use==conditionsLevel[1]))]
        dataCon2 <- count_norm_use[,c(which(conditions_use==conditionsLevel[2]))]
        foldChanges <- log2(rowMeans(dataCon2)/rowMeans(dataCon1))

        ## Output results based on FDR threshold
        outRst<-data.frame(log2foldChange=foldChanges, pValues=pvalues, FDR=fdr)
        rownames(outRst)=rownames(count_norm_use)
        outRst=na.omit(outRst)

        ## class1代表分子
        outRst$class1 <- class2
        outRst$class2 <- class1
        outRst$meanTMM_class1 <- rowMeans(dataCon2)
        outRst$meanTMM_class2 <- rowMeans(dataCon1)
        outRst$gene_id <- rownames(outRst)

        result_diff <- rbind( result_diff , outRst)
    }
}

##########################################################################################

colnames(result_diff) <- c("log2FoldChange" , "pvalue" , "padj" , "class1" , "class2" , "meanTMM_class1" , "meanTMM_class2" , "gene_id")

image_name <- paste0( out_path , "/DiffGene.tsv" )
write.table( result_diff , image_name , row.names = F , sep = "\t" , quote = F )
